\section{Diffusion Tensor Imaging}

Diffusion tensor Imaging (DTI) is a technique to measure the anisotropic diffusion properties of biological tissues within the sample. This allows us to noninvasively infer the structure of the underlying tissue. Diffusion properties allow the classification of different types of tissues and can be used for tissue segmentation and detecting tissue orientations. The principal eigenvector of the diffusion tensor is known to align with fiber tracts in the brain \cite{pierpaoli96} and in the heart \cite{scollan98}. Heart fibers reconstructed from cardiac DTI are shown in Figure \ref{fibers-dti}.

Diffusion tensors describe the diffusion properties of water molecules. In tissues diffusion properties are dictated by the cell structure of the tissue. Since cell membranes are selectively permeable, water molecules can move easily within a cell, but their diffusion across the membrane is limited. Thus diffusion properties of the tissue reflect the shape and orientation of the cells. For the specific case of elongated cells like the cardiac muscles, the diffusion will be maximum along the primary axis of the muscle, which also happens to be the direction along which maximum strain is developed.

\begin{figure}
\begin{center}
\includegraphics[width=0.4\textwidth]{images/fibers3}
\caption{Heart fiber orientation in the human heart, obtained from diffusion tensor imaging}
\label{fibers-dti}
\end{center}  
\end{figure} 

Diffusion is measured through a diffusion coefficient, which is represented as a symmetric second order tensor:

\begin{equation}
\label{eq:dt}
{\bf D} = \left( 				
\begin{array}{ccc}
D_{xx} & D_{xy} & D_{xz} \\
D_{yx} & D_{yy} & D_{yz} \\	
D_{zx} & D_{zy} & D_{zz}
\end{array}
          \right)
\end{equation}
The 6 independent values of the tensor elements vary continuously with the spatial location in the tissue.

Eigenvalues $\lambda_i$ and eigenvectors ${\bf e}_i$ of the diffusion tensor (\ref{eq:dt}) can be found
as a solution to the eigenvalue problem:
\[
{\bf De}_i = \lambda_i{\bf e}_i
\]
Since the tensor is symmetric, its eigenvalues are always real numbers, and the eigenvectors are orthogonal and form a basis. Geometrically, a diffusion tensor can be thought of as an ellipsoid with its three axes oriented along these eigenvectors, with the three semiaxis lengths proportional to the square root of the eigenvalues of the tensor.

% Using the ellipsoidal interpretation, one can classify the diffusion properties of a tissue according to the shape of the ellipsoids, with extended ellipsoids corresponding to regions with strong linear diffusion (long, thin cells), flat ellipsoids to planar diffusion, and spherical ellipsoids to regions of isotropic media (such as fluid filled regions like the ventricles). The quantitative classification can be done through the coefficients $c_l$, $c_p$, $c_s$ corresponding to linear, planar and spherical diffusion.
%\begin{eqnarray*}
%c_l &=& \frac{\lambda_1 - \lambda_2}{\lambda_1 + \lambda_2 + \lambda_3} \\
%c_l &=& \frac{2(\lambda_2 - \lambda_3)}{\lambda_1 + \lambda_2 + \lambda_3} \\
%c_l &=& \frac{3\lambda_1}{\lambda_1 + \lambda_2 + \lambda_3}
%\end{eqnarray*}

The heart fibers have strong linear diffusion %($\lambda_1 >> \lambda_2 \approx \lambda_3)$ 
and are oriented along the principal eigenvector, ${\bf e}_1$ \cite{tseng99, scollan98}. Therefore if we need fiber orientations for a specific subject, then computing the principal direction of the diffusion tensor is sufficient. However, {\em in vivo} acquisition of cardiac DTI is not possible using current scanners and acquisition protocols. Consequently, we use an alternate strategy and map the diffusion tensors from a template image onto the subject. The procedure for warping DTIs from a template onto a subject is described in the following section.

\section {Warping Diffusion Tensors from template to subjects}
\label{sec:warpDTI}

We have MR images for both the subject and the template. The diffusion tensor data is available only for the template. We use a very high dimensional elastic registration technique \cite{hammer} to estimate the deformation that warps the template to the subject space. We use this deformation field to map the fibers from the template to the subject. It is more complicated to warp tensor fields than it is to warp scalar images. This is because the tensor must be reoriented on each image voxel, in addition to a voxel displacement that is implied by the deformation field. This is achieved by finding the rotational component of the deformation field. We test the effectiveness of the tensor remapping algorithm by comparing the mapped tensors with the ground truth diffusion tensors for 19 canine datasets\footnote{CCBM, Johns Hopkins University}. The method of computing the transformation between the two geometries and the tensor reorientation algorithm are now described.

\subsection{Deformable Image Registration}

Image warping for deformable registration has received a great deal of attention during the past decade \cite{Zitova03}. In the present work we used a very-high-dimensional elastic transformation procedure in 3D volume space, referred to as the hierarchical attribute matching mechanism for elastic registration (HAMMER) method, which is determined from T1-weighted images and applied on the coregistered DT image of the template. This approach uses image attributes to determine point correspondences between an individual image and a template, which resides in the stereotaxic space and is the subject for which we have the diffusion tensors. A hierarchical sequence of piece-wise smooth transformations is then determined, so that the attributes of the warped images are as similar as possible to the attributes of the target. Relatively fewer, more stable attributes are used in the initial stages of this procedure, which helps avoid local minima, a known problem in high-dimensional transformations. The details of this algorithm can be found in \cite{hammer}.


\subsection{Tensor Reorientation}
It is a simple matter to warp a scalar image by a known spatial transformation. The image value from a particular voxel is transferred, via the displacement field of the spatial transformation, to a voxel in the target image. Typically, some sort of interpolation must also be applied. However, a more complex procedure is required to warp tensor fields, especially when the tensor estimates are noisy. We use an approach similar to that proposed by Xu et al \cite{xu03}.

If we know the direction, ${\bf v}$, of the fiber on voxel with coordinates ${\bf x}$, we can readily find the rotated version, ${\bf v}'$, of ${\bf v}$, according to the warping transformation. If ${\bf R}$ is the matrix that rotates ${\bf v}$ to ${\bf v'}$, then ${\bf R}$ should be applied to the respective tensor measurement. However, in practice we do not know ${\bf v}$. In fact, this is precisely what we would like to estimate. We only have a noisy orientation of $v$, which is the principal direction (PD) of the corresponding tensor measurement. One could use that PD in place of ${\bf v}$, as proposed in \cite{alex01}. However, that makes the approach vulnerable to noise, since the PD is only a noisy observation, and could be quite different from the true underlying fiber orientation. 

Assuming that we know the probability distribution function (PDF), $f({\bf v})$, of the fiber direction ${\bf v}$, we can find the rotation matrix, $\tilde{\bf R}$ which minimizes the expected value of $\|{\bf v' -Rv} \|^2$ over all orthonormal matrices ${\bf R}$:

\begin{eqnarray*}
\tilde{\bf R} &=& \mathop{\arg \min}_{\bf R} E\{\|{\bf v' -Rv} \|^2\}\\
&=& \mathop{\arg \min}_{\bf R} \int_{\bf v} pdf({\bf v}) \|{\bf v' -Rv} \|^2 d{\bf v}
\end{eqnarray*}

This problem can be solved by the Procrustean estimation \cite{golub83}, if a number of random samples, ${\bf v}$, are drawn from the PDF, and their respective rotated versions, ${\bf v'}$, are found by the rotation that the warping field applies to ${\bf v}$. If we arrange these vectors ${\bf v'}$ and ${\bf v}$ to form the columns of the matrices ${\bf A}$ and ${\bf B}$, respectively, then $\tilde{\bf R}$ is found by minimizing:
\[
	\left\|{\bf A} - {\bf R}\cdot{\bf B}\right\|_2^2 = \left\|{\bf A}\right\|_2^2 + \left\|{\bf B}\right\|_2^2 + 2\sum_i\sigma_i({\bf A}\cdot{\bf B}^T)
\]

where $\sigma_i(\bM)$ are the singular values of matrix ${\bf M}$. $\tilde{\bf R}$ can be determined via a singular value decomposition of ${\bf A \cdot B}$.

\begin{eqnarray*}
{\bf A}\cdot{\bf B}^T &=& {\bf V}\cdot\Sigma\cdot {\bf W}^T \\
\tilde{\bf R} &=& {\bf V}\cdot{\bf W}^T 
\end{eqnarray*}

More details on the algorithm and the estimation of the PDF can be found in \cite{xu03}.
